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Abstract. 

The anisotropic phase-field crystal model recently proposed and used by Prieler 
et al. [J. Phys.: Condens. Matter 21, 464110 (2009)] is derived from microscopic 
density functional theory for anisotropic particles with fixed orientation. Further 
its morphology diagram is explored. In particular we investigated the influence of 
anisotropy and undercooling on the process of nucleation and microstructure formation 
from atomic to the microscale. To that end numerical simulations were performed 
varying those dimensionless parameters which represent anisotropy and undercooling 
in our anisotropic phase-field crystal (APFC) model. The results from these numerical 
simulations are summarized in terms of a morphology diagram of the stable state 
phase. These stable phases are also investigated with respect to their kinetics and 
characteristic morphological features. 



PACS numbers: 81.10.Aj, 64.70.dm, 82.70.Dd 
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1. INTRODUCTION 

Early in the 90th Swift and Hohenberg formulated an amplitude approach to describe 
systems, where the stable states are periodic, as e.g. the case for Rayleigh-Benard 
convection [1]. More recently this idea has been taken up by the materials science 
community to model crystals at the atomic scale. Elder et al. proposed a functional for 
a scalar dimensionless field of form 

P = J{1^ M + ^'f ^0^) dr , (1) 

V 

with two phenomenological parameters go and e and a corresponding dynamical equation 




for this purpose [2]. Since its introduction, this phase-field crystal (PFC) method 
[21 [31 m El E] has emerged as a computationally efficient alternative to molecular 
dynamics (MD) simulations for problems where the atomic and the continuum scale are 
tightly coupled. The reason is that it operates for atomic length scales and diffusive time 
scales. Thus for a simple application such as diffusion in gold or copper it runs 10^-10^ 
times faster than the corresponding MD calculation [7]. In that sense it provides from 
point of view of multiscale materials modeling an interesting link between the traditional 
phase-field method and MD. Moreover, a connection between classical density functional 
theory of freezing and phase- field crystal modeling could be identified in [4l[8]. Thereby a 
second theoretical foundation besides the Swift-Hohenberg amphtude equation approach 
could be established. Essentially it motivates the application of PFC models also for 
spatially non-uniform non-periodic states. 

Recently the phase-field crystal method has been applied to a variety of different 
growth phenomena. One of its interesting features is that other than the phase-field 
method, in which elasticity explicitly needs to be integrated in the functional to be 
taken into account [9j|, it includes elastic effects inherently. Thus it allows to simulate 
for example features of crack propagation [1] and plasticity [31 HO] from the atomic 
to the micro-scale. To model the elastic behavior of different kinds of materials, the 
parameters of the phase-field crystal model equation can be adjusted to match the elastic 
moduli of a given experimental system. However, in its most simplistic form, in which 
it is a reformulation of the Swift-Hohenberg equation [1] with a conserved dynamics as 
introduced by Elder et al. [21 13] the Poisson ratios which can be modeled, are restricted 
to 1/3 (in the one mode approximation). Moreover, since in the simplistic PFC model 
the material is defined by only three parameters, it is restricted with respect to the 
crystal lattice structures which it can describe as well. These are triangular symmetries 
in two dimensions and BCC symmetry in three dimensions [HI |12]. Another crystal 
symmetry applying to protein crystals in a membrane could be obtained by including 
higher order correlation functions jlS]. Moreover, liquid crystals have been simulated by 
combining the original phase- field crystal equation with an orientational field [T^fT6l fT5]. 
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In [T3| we followed the above direction to extend the phase-field to apply to a larger 
class of condensed matter systems taking a different route: We derived a generalized 
PFC model for isotropic as well as anisotropic crystal lattice systems of arbitrary Poisson 
ratios as well as condensed matter systems built-up from non-spherical units such as for 
example anisotropic colloids and liquid crystals. To this end we extended the simplest 
PFC model proposed by Elder [2], to a conservative, anisotropic Langevin equation and 
applied it to the heterogeneous nucleation of colloids at a wall. 

Whereas our previous work [TTj was devoted to a mere introduction of our model 
and its application to heterogenous nucleation at a wall, here we show for the first time 
how its parameters can be derived from dynamical density functional theory (DDFT) 
(see section |2]). Further we report for the first time in detail on its morphology diagram. 
To do so, we proceed as follows: First we give a thorough derivation of our anisotropic 
phase-field model based on the DDFT in section [21 We then - after a brief summary 
of the model in section |2] - study the infiuence of anisotropy and undercooling on the 
morphology of the final states. We analyze these morphologies and summarize our 
results in terms of a morphology diagram in section HJ Finally we conclude with a 
summary and an outlook of our study in section [5l 



2. ANISOTROPIC PHASE-FIELD CRYSTAL MODEL 



2.1. The model 

By now phase-field crystal (PFC) modeling is widely used to predict crystal nucleation 
growth and to model microstructural pattern formation during different physical 
phenomena such as solidification. As usual, the PFC model used in this studies is based 
on a free energy fuctional F[(j){r, t)] of phase-field (j){r, t) and a dynamical equation which 
represents the time evolution of the phase-field. In this PFC model, periodic nature of 
a crystal lattice is incorporated by using a free energy functional which is minimized by 
periodic density field. The equation of motion used in this model was introduced in [2] 
for the case of simplest phase-field crystal (SPFC) model, and is given by 

= A[{(go^ + A)^-e}0 + 0^], (3) 

here, Qq and e are constants. In order to simplify the model, a dimensionless parameter 
T is introduced which is defined as 

r = -{ql-e). (4) 



In our model r represents undercooling in the same manner as r is defined in jl8] . 
therefore, r can be written as 

r oc AT. (5) 

The anisotropic version of phase-field crystal model originally introduced in [19] is used 
in this study. This PFC model is basically an extension of the SPFC model which 
is derived in [20j. The APFC model is capable of simulating isotropic and anisotropic 
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crystal lattice system of any arbitrary poisson ratio as well as condensed matter systems 
such as colloids and liquid crystals. The free energy functional used in this model is 
given by 



T + a. 



+ h 



ijkl 



dxidxjdxkdxi 



+ \ccP^)dr, (6) 



where dij IS cL symmetric matrix and bijti is a fourth rank tensor with the symmetry of 
an elastic tensor: i ^ j, k -v^ I, {i,j) ^ {k, I). From the free energy functional defined in 
Equation ([6]), the corresponding Langevin differential equation of motion for anisotropic 
lattice system can be written as follows: 

T + aijT: — H bijki- — - — - — — (p + cq 



P 



dt 



A 



dxidxj 



' dxidxjdxkdxi- 



(7) 



2.2. Derivation of the anisotropic phase- field crystal model from dynamical density 
functional theory 

The coefficients occuring in the anisotropic phase-field crystal model proposed in Ref. 
pT] can be derived from microscopic density functional theory [211 1221 [231 121] • Here we 
follow a similar line as proposed recently by van Teeffelen et al |H] for radially symmetric 
interactions. We generalize this route here to anisotropic interactions. 

We assume that the anisotropic colloids are completely aligned in space. Cartesian 
coordinates r = {xi,X2, ...,Xd) will be used in the following, d denoting the spatial 
dimension. The interaction pair potential between two aligned particles is M(r)[25]. The 
latter function is anisotropic, in general, i.e. it does not only depend on |r|. Other 
examples for these anisotropic interactions with fixed orientations are oriented hard 
spherocylinders [26] and charged rods [271 EH], anisotropic Gaussian potentials [29] , 
board- like colloidal particles [30], colloidal molecules [31], as well as patchy colloids [32] 
and proteins [HH [34] • Henceforth inversion symmetry is assumed 

u{-r) = u{r) (8) 

Dynamical density functional theory for anisotropic situations [35] is now generalized 
from the isotropic case as follows. The dynamical evolution of the time-dependent one- 
particle density field p(r, t) is: 

,5F[piv,t)] 



p{r,t) = {kBTy'V ■ 



Vp{r,t)V- 



(9) 



<5p(r,t) J 

Here fc^T is the thermal energy, and V = {d/dxi, d/dx2, d/dxd) is the d-dimensional 
gradient. V = diag{Di, D2, Dd) denotes the diagonalized diffusion tensor with the 
anisotropic short-time translational diffusivities of the anisotropic particle. For a given 
(hydrodynamic) shape of the particle, explicit expressions for Di are available [511 ISZ] ■ 
Furthermore, in Equation Q, F [p(r, t)] is the equilibrium density functional which can 
be split as 

F [p(r)] = Fid [p(r)] + F,^ [p(r)] + F,^, [p(r)] . (10) 
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Fid [p(r)] = ksT I drp(r) {in [p(r)A'^] - l} , (11) 

with A denoting the thermal de Broglie wavelength. The external part involves an 
external one-body potential V{r,t) and is given by 

Fe.t[p(r)] = j drp{r)V{r,t) . (12) 

Finally, the excess part Fex[p(r)], embodies the nontrivial correlations between the 
particles and must be further approximated. Henceforth we assume small deviations 
of the inhomogeneous density profile around a homogeneous reference density p. In 
this limit, the leading approximation for [p(r)] is given by the Ramakrishnan and 
Yussouff [38] expression: 

Fox [p(r)] ^ FM -^11 drdr'Ap(r)Ap(r')4^)(r - r'; p) . (13) 

(2) 

where Cq (r — r'; p) is the anisotropic direct correlation function of the fluid at density p 
which possesses the same symmetry as the underlying pair potential u{r). In particular, 
it is inversion-symmetric 

4'^(-r,p) = -4'^(r,p) (14) 

Moreover, Ap(r) = p(r) — p. In Fourier space Equation (fT3|) reads 

^ex[p(r)] = FUp) - I dkAp(k)Ap(-k)4^)(k,p) (15) 

with ~ denoting a Fourier transform. We now expand the direct correlation function 
Cq (k, p) in terms of k around k = 0. (Alternatively fitting procedures can be used, e.g. 
around the first peak of CQ^''(k, p).) This leads to the Taylor expansion in Fourier space 

d d 

CQ\k, p) = Co + ^ aijkikj + ^ hjkikikjkkh + ... . (16) 

i,j=l i,j,k,l=l 

corresponding to a gradient expansion in real-space. Inversion symmetry Equation (1141) 
enforces all odd orders to vanish. Possible additional symmetries in the shape of the 
particles will lead to corresponding restrictions on the tensorial coefficients aij and bijki 
as discussed below. 

Inserting this expansion into Equation ([9]), one gets 



p(r, t) = V ■ Wp(r, t) + V-VV ^{kBT)-'V{r, t) 



d d „4 



JJ = 1 i,j,k,l=l ■' 

If one further uses the constant mobility approximation, p(r, t) = p in front of the 
functional derivative in Equation (jH]) and if one approximates 

F,a[pir)]^kBTpJ dr|i0(r,t)2-i0(r,t)3 + l0(r,t)4- const.} (18) 
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with (p{r,t) = Ap{r,t)/p, one arrives at: 
0(r,t) =pV-PV 



r, t) - ^0(r, tf + i0(r, tf + (fc^T)- V(r, t) 



i,j=l i,j,k,l=l ■' 



(19) 



This exactly reduces to the anisotropic phase-field model of Ref. [IT] for the special case 
d = 2, V = Dqoo, and a neglected cubic term in the ideal gas functional expansion in 
Equation (ITSll . As a remark the latter was retained in other variants of the PFC model 



Concluding this section, the anisotropic phase-field crystal model as used in flT] can 
be derived and justified from dynamical density functional theory. The derivation points, 
however, to more realistic approximations for anisotropic diffusivities. Furthermore, if 
Equation (ITSil is used, some approximations can be avoided but these were not found 
to change the results significantly for spherical interactions [8]. 

2.3. Phenomenological symmetry considerations 

We finally present phenomenological symmetry arguments for the expansion coefficients 
Qij and hijki of the anisotropic PFC model. First we assume that the orientation of 
the fixed particles is set by a single unit vector E only which is invariant under space 
inversion (r — )■ — r). This is the case for d = 2 and for rotationally symmetric particles 
in d = 3. Then, any gradient term in the scalar free energy functional must involve 
an even number of gradients due to space inversion symmetry. Rotational symmetry of 
space then requires that only combinations of ^ ■ V and V ■ V are nonvanishing in the 
functional. Therefore the only possibility for physically relevant gradient terms is 

X:a,-^ = Ai(^.vr + A,A (20) 

i,j=l * ^ 

and 

i,j,k,l=l 

where Ai, A2, A3, A4, and A5 are scalar prefactors. This reduces the number of 
independent degrees of freedom in aij and hijki down to 5. 

In case there are different fixed vectors, say E and 5, there are corresspondingly 
more terms allowing for more freedom in a^j and hij^i- This is realized, e.g., for biaxial 
colloidal particles in two crossed external fields along E and B. 



3. SIMULATION PARAMETERS 



As initial condition, a square domain is defined with a sphere in the center to initialize 
the nucleus. Periodic boundary conditions in all directions in the square box are used. 
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The values for the parameters are the same as defined in [T7] , namely p = I (which sets 
the time scale), an = = 2, buu = 62222 = 61122 = 1, 61212 = and c = 1, respectively 
[H]. However, a typical set of values for r and s is used for each simulation since our 
basic objective is to study the dependence of the stable state phase on these parameters. 
A simple explicit numerical scheme is used to obtain a reasonably well approximated 
solution. A forward Euler scheme is used for the time derivative with a sufficiently small 
time step of At =0.00075 to ensure the stability of the scheme. The Laplace operator 
is approximated by using second order difference scheme given by 



= (^(pi+ij + (pi^ij + (pij+i + 0ij_i - 40i,,)/(Ax)2. (22) 

For the following simulations Ax is chosen as 11/4. Convergence of our results was 
ensured via convergence studies. The morphologies depicted in figure [H figure [2] and 
figure m show simulations of 256 times 256 numerical grid units. However, morphologies 
in figure [5] show simulation of 128 times 256 numerical grid units. 



4. RESULTS AND DISCUSSION 



In this section we present the simulation results obtained from our studies. These 
simulation results demonstrated the following issues concerning nucleation and 
successive microstructure formation: 

(i) The effect of undercooling on crystal growth. 

(ii) The dependence of anisotropy and undercooling on stable state phase. 

(iii) The effect of anisotropy and undercooling on distance between the neighbouring 
stripes. 



4.1. Anisotropic effects 

In order to quantify the anisotropy of the material at the atomic scale, a dimensionless 
parameter s is introduced and is given by 

^1112 /^o\ 

s = -- . (23) 

61111 

The effect of this dimensionless parameter is studied by performing numerical 
simulations with s = and s = 0.125. The initial and boundary conditions as well 
as the values for all other simulation parameters used in these simulations as same as 
given in section [3l The results obtained for both cases after 30,000 time steps are shown 
in figure [H 



4-2. Undercooling effects 

When a liquid is supercooled just below the melting temperature the crystal starts 
growing and the crystal growth is directly related to the undercooling. Depending on 
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a) 




b) 



Figure 1. Simulation results of crystal growth calculated for r = —3/4 and (a) s = 
(isotropic case) and (b) s — 0.125 (anisotropic case). All other parameters are chosen 
as given in section [31 As expected the isotropic case shows a symmetric morphology. 

the formal undercooling, which quantifies the distance from the phase-equilibrium line 
in the phase diagram, in the system the final state will have different morphologies. 
These are categorised and analysed here in detail to get an overview of the state phase 
of the APFC model given by the variables s and r. 

In this section the APFC model is used to examine the rate of crystal growth from 
a supercooled liquid state. As explained above, r represents the undercooling in our 
model. A number of simulations with different values of r are performed for a specific 
s value. We used the same initial condition i.e. single solid nucleus (nucleation site) for 
all simulations. The results showed that the rate of crystal growth increases with an 
increase in the value of r. As an example the simulation results with s = and for two 
typical values of r i.e. with r = —0.25 and r = —0.8 respectively at 40,000 time steps 
are depicted in figure El 




Figure 2. Simulation results of crystal growth after 40,000 time steps for s = and 
(a) T = -0.8 and (b) r = -0.25 

4-3. Heterogeneous nucleation and crystal growth 

In this section we studied the dependence of a stable state phase on anisotropy and 
undercooling. More specifically, we show how a stable state can be composed of different 
phases such as triangular phase, stripe phase and co-existance of stripes and triangular 
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phase, depending on the values of s and r i.e. anisotropy and undercoohng respectively. 
To study these patterns we performed a number of simulations with different values of 
s and r in each simulation. However, the other parameters for all these simulations are 
same as described in the previous section. The results in the form of a diagram of stable 
state phase are demonstrated in figure [3l 
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Figure 3. Simulation results of stable phases for same initial conditions and different 
values of s and t 

These simulation results demonstrate that the stable state phase always consists of 
stripes if s > 0.25 irrespective of the r value. The co-existence of stripes and triangular 
phases is found only in case of s = 0.125 and r > —0.5, while for other values of s and 
r, the stable state consist of a triangular phase. Typical picture of the triangular phase, 
stripes, and a co-existence of stripes and triangular phase are given in figure |H 

The precise shape of the grain boundary between the two solid phases depends on 
the system size due to |12] . The two-phase coexistence as such, however, is independent 
of the size of the system. 

4-. 4- Distance between the neighbouring stripes 

To analyse the stripe morphology further, we studied the effect of anisotropy and 
undercooling on the distance between the neighbouring stripes. As discussed in the 
previous section, the stable state consists of stripes when s > 0.25 irrespective of the r 
value, as shown in figure |3l It is observed that the stripe phase obtained for different 
values of s are different from each other in terms of the spacing between the neighbouring 
stripes in the stripe phases. However, simulations with different values of r result in 
similar stripe phases. Our results reveal that the spacing between neighbouring stripes 
decreases with an increase in the value of s. Thus the stripe phase obtained with s = 1 
are much finer compared to the ones obtained with s = 0.25. Further, to investigate the 
effect of the undercooling r on the stripe phase, we performed several simulations by 
fixing a specific value for s and varying r. The results show that for a specific value of s, 
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Figure 4. Typical stable states phases from simulations performed for (a) s = 0.125 
and r = —0.25, which result in a co-existence of stripes and triangular phase, (b) s = 
and T = —0.75 which result in the triangular phase , and (c) s — 0.75 and t = —0.25, 
which results in stripes 

similar stripe morphologies are obtained with different value of r. Figure [5] demonstrates 
this interesting finding for different values of s and r. 

— I 1 1 1 1 — 

■ ■■iii 




0.25 0.5 0.75 1 

S 

Figure 5. Stripe phases for different values of s and r. Resulting stripe phase for (a) 
s = 0.25, (b) s = 0.625 and (c) s = 1. 

From the above discussion, we can conclude that the distance between neighbouring 
stripes decreases as we increase the value of s. However, the undercooling has no 
significant effect on the stable stripe phase, i.e. the stripe phases for different values of 
r are similar. 
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4-5. Analysis of the anisotropy introduced by s 

In this section we discuss the influence of the anisotropic factor s on the final triangular 
phase obtained with < s < 0.125 and r = —0.75. More specifically, we analysed 
the form of triangles in each of the final triangular phase. The form of a triangle is 
determined in terms of three internal angles and ratio of length of the longest and 
shortest sides of the triangle. It is observed that the triangular phase obtained in case 
of s = i.e. without any anisotropy, consist of triangles with same three angles of 
60° each. However, the triangular phase obtained with non-zero values of s contains 
triangles with dissimilar sides. The details of the angles calculated for each case are 
shown in figure El 
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Figure 6. Form of triangles in the final triangular phase obtained for different 
anisotropies s. 

These results demonstrate that for isotropic case i.e. s = 0, the final triangular 
phase consist of equilateral triangles. However, for anisotropic case i.e. non zero values 
of s, the final triangular phase contain scalene triangles i.e. no two sides are similar. 
The ratio of length of the longest and shortest sides of the triangles, calculated for each 
case is presented in figure [71 This underlines the capability of our APFC model to give 
rise to truly anisotropic morphologies. 

5. SUMMARY AND OUTLOOK 

In this article we presented a DDFT based derivation of the APFC model proposed by 
two of the authors (D. L. and H. E.) in [T7] previously. Further we investigated the 
state phase of this model given by variation of r and s to demonstrate its capacity to 
model structures beyond those captured by the SPFC equations originally introduced 
by Elder et al. [2]. 

In particular we studied the influence of anisotropy and undercooling on final states 
using numerical techniques to minimize the free energy functional in our model. More 
specifically, a number of numerical simulations are performed by using different sets 
of values for our dimensionless parameters s and r which represent anisotropy and 
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undercooling, respectively. The results obtained from these numerical simulations are 
analysed. Our studies reveal that: 

(i) the rate of crystal growth increases with increase of r i.e. undercooling; 

(ii) the stable state phase consists of a stripe phase if s > 0.25 irrespective of r. However, 
the stable state is a co-existence of stripes and triangular phase when s = 0.125 and 
T > —0.5, while for other values of s and r, the stable state consists of a triangular 
phase; 

(iii) for s > 0.25, the undercooling r has no effect on the resulting stripe phases, however, 
the distance between the neighbouring stripes decreases with an increase of s. 

(iv) Triangular crystals with cells that are neither equilateral nor rhombic are possible 
for our anisotropic model. 

In the future we plan to extend the approach to reactive systems [H] to simulate 
morpholgenesis in such systems from the atomic to microscale. 
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